#include "stdafx.h"
#include "NetAlign.h"
#include <omp.h>

#define			FirstIter			0
#define			SecondIter		9
#define			FNodeNum			4000
#define			GNodeNum			4000

TIntSet FSet, GSet;
TIntH FMatchH, GMatchH;
TIntFltH FPageRankH, GPageRankH;

TIntH FNodeIdH, GNodeIdH;
double CurR[FNodeNum][GNodeNum], NextR[FNodeNum][GNodeNum];
int DegF[FNodeNum], DegG[GNodeNum];
int LenF[FNodeNum], LenG[GNodeNum];
int NbF[FNodeNum][FNodeNum], NbG[GNodeNum][GNodeNum];
short int NbSF[FNodeNum][FNodeNum], NbSG[GNodeNum][GNodeNum];
double C[FNodeNum][GNodeNum];

void ComputeFDegNb(PNGraph& F, TStr Pref) {
	for (TNGraph::TNodeI NI = F->BegNI(); NI < F->EndNI(); NI++) {
		int Id = NI.GetId();
		FNodeIdH.AddDatId(Id);
	}

	memset(LenF, 0, sizeof(LenF));
	memset(DegF, 0, sizeof(DegF));
	for (TNGraph::TNodeI NI = F->BegNI(); NI < F->EndNI(); NI++) {
		int SrcId = NI.GetId(), SrcN = FNodeIdH.GetDat(SrcId).Val;
		TIntSet Visited;
		for (int e = 0; e < NI.GetOutDeg(); e++) {
			int DstId = NI.GetOutNId(e), DstN = FNodeIdH.GetDat(DstId).Val;
			Visited.AddKey(DstN);
			DegF[SrcN] += 10;
		}
		for (int e = 0; e < NI.GetInDeg(); e++) {
			int DstId = NI.GetInNId(e), DstN = FNodeIdH.GetDat(DstId).Val;
			if (Visited.IsKey(DstN)) continue;
			DegF[SrcN] += 1;
		}
	}

	for (TNGraph::TNodeI NI = F->BegNI(); NI < F->EndNI(); NI++) {
		int SrcId = NI.GetId(), SrcN = FNodeIdH.GetDat(SrcId).Val;
		TIntSet Visited;
		for (int e = 0; e < NI.GetInDeg(); e++) {
			int DstId = NI.GetInNId(e), DstN = FNodeIdH.GetDat(DstId).Val;
			Visited.AddKey(DstN);
			NbF[SrcN][LenF[SrcN]] = DstN;
			NbSF[SrcN][LenF[SrcN]] = 10;
			LenF[SrcN] += 1;
		}

		for (int e = 0; e < NI.GetOutDeg(); e++) {
			int DstId = NI.GetOutNId(e), DstN = FNodeIdH.GetDat(DstId).Val;
			if (Visited.IsKey(DstN)) continue;
			NbF[SrcN][LenF[SrcN]] = DstN;
			NbSF[SrcN][LenF[SrcN]] = 1;
			LenF[SrcN] += 1;
		}
	}
	FILE* FLog = fopen(Pref.CStr(), "w");
	for (int i = 0; i < FNodeNum; i++) {
		fprintf(FLog, "Node %d : Deg = %f\n", i, DegF[i]);
	}

}

void ComputeGDegNb(PNGraph& G, TStr Pref) {
	for (TNGraph::TNodeI NI = G->BegNI(); NI < G->EndNI(); NI++) {
		int Id = NI.GetId();
		GNodeIdH.AddDatId(Id);
	}

	memset(LenG, 0, sizeof(LenG));
	memset(DegG, 0, sizeof(DegG));
	for (TNGraph::TNodeI NI = G->BegNI(); NI < G->EndNI(); NI++) {
		int SrcId = NI.GetId(), SrcN = GNodeIdH.GetDat(SrcId).Val;
		TIntSet Visited;
		for (int e = 0; e < NI.GetOutDeg(); e++) {
			int DstId = NI.GetOutNId(e), DstN = GNodeIdH.GetDat(DstId).Val;
			Visited.AddKey(DstN);
			DegG[SrcN] += 10;
		}
		for (int e = 0; e < NI.GetInDeg(); e++) {
			int DstId = NI.GetInNId(e), DstN = GNodeIdH.GetDat(DstId).Val;
			if (Visited.IsKey(DstN)) continue;
			DegG[SrcN] += 1;
		}
	}

	for (TNGraph::TNodeI NI = G->BegNI(); NI < G->EndNI(); NI++) {
		int SrcId = NI.GetId(), SrcN = GNodeIdH.GetDat(SrcId).Val;
		TIntSet Visited;
		for (int e = 0; e < NI.GetInDeg(); e++) {
			int DstId = NI.GetInNId(e), DstN = GNodeIdH.GetDat(DstId).Val;
			Visited.AddKey(DstN);
			NbG[SrcN][LenG[SrcN]] = DstN;
			NbSG[SrcN][LenG[SrcN]] = 10;
			LenG[SrcN] += 1;
		}
		for (int e = 0; e < NI.GetOutDeg(); e++) {
			int DstId = NI.GetOutNId(e), DstN = GNodeIdH.GetDat(DstId).Val;
			if (Visited.IsKey(DstN)) continue;
			NbG[SrcN][LenG[SrcN]] = DstN;
			NbSG[SrcN][LenG[SrcN]] = 1;
			LenG[SrcN] += 1;
		}
	}
	FILE* FLog = fopen(Pref.CStr(), "w");
	for (int i = 0; i < GNodeNum; i++) {
		fprintf(FLog, "Node %d : Deg = %f\n", i, DegG[i]);
	}
}

void InitializeC(PNGraph& M) {
	memset(C, 0, sizeof(C));
	TIntSet FTrueSet, GTrueSet;
	for (TNGraph::TEdgeI EI = M->BegEI(); EI < M->EndEI(); EI++) {
		int FId = EI.GetSrcNId(), FN = FNodeIdH.GetDat(FId).Val;
		int GId = EI.GetDstNId(), GN = GNodeIdH.GetDat(GId).Val;
		C[FN][GN] = 1;
		FTrueSet.AddKey(FN);
		GTrueSet.AddKey(GN);
	}

	double Z = GNodeNum - FTrueSet.Len();
	for (int i = 0; i < FNodeNum; i++) {
		if (FTrueSet.IsKey(i)) continue;
		for (int j = 0; j < GNodeNum; j++) {
			if (GTrueSet.IsKey(j)) continue;
			C[i][j] = 1.0 / Z;
		}
	}
}

void InitializeR() {
	for (int i = 0; i < FNodeNum; i++)
		for (int j = 0; j < GNodeNum; j++)
			CurR[i][j] = 1.0 / GNodeNum / DegF[i] / DegG[j];
}

double ComputeRelativeRDif() {
	double Dif = 0, Sum = 0;;
	for (int i = 0; i < FNodeNum; i++)
		for (int j = 0; j < GNodeNum; j++) {
			double dist = (CurR[i][j] * DegF[i] * DegG[j] - NextR[i][j]);
			Dif += dist * dist;
			Sum += NextR[i][j] * NextR[i][j];
		}
	Dif = sqrt(Dif / FNodeNum);
	Sum = sqrt(Sum / FNodeNum);
	printf("Dif = %lf, Sum = %lf, Error Percentage = %lf%%\n", Dif, Sum, 100*(Dif/Sum));
	return (Dif/Sum);
}

void ComputeC(PNGraph& M, double thresh, TStr Pref) {
	memset(C, 0, sizeof(C));
	TIntSet FTrueSet, GTrueSet;
	for (TNGraph::TEdgeI EI = M->BegEI(); EI < M->EndEI(); EI++) {
		int FId = EI.GetSrcNId(), FN = FNodeIdH.GetDat(FId).Val;
		int GId = EI.GetDstNId(), GN = GNodeIdH.GetDat(GId).Val;
		C[FN][GN] = 1;
		FTrueSet.AddKey(FN);
		GTrueSet.AddKey(GN);
	}

	TIntFltH RowMinThH, ColMinThH, RowThH, ColThH;
	// Select top candidates in each row
	for (int i = 0; i < FNodeNum; i++) {
		if (FTrueSet.IsKey(i)) continue;
		TFltV ValV;
		for (int j = 0; j < GNodeNum; j++) {
			if (GTrueSet.IsKey(j)) continue;
			C[i][j] = CurR[i][j] * DegF[i] * DegG[j];
			ValV.Add(C[i][j]);
		}
		ValV.Sort(false);
		double ThreshVal = ValV[int(FNodeNum * thresh)].Val;
		RowMinThH.AddDat(i) = ThreshVal;
		RowThH.AddDat(i) = ValV[2];
	}

	for (int j = 0; j < GNodeNum; j++) {
		if (GTrueSet.IsKey(j)) continue;
		TFltV ValV;
		for (int i = 0; i < FNodeNum; i++) {
			if (FTrueSet.IsKey(i)) continue;
			ValV.Add(C[i][j]);
		}
		ValV.Sort(false);
		double ThreshVal = ValV[int(FNodeNum * thresh)].Val;
		ColMinThH.AddDat(j) = ThreshVal;
		ColThH.AddDat(j) = ValV[2];
	}

	// Filtering
	for (int i = 0; i < FNodeNum; i++) {
		if (FTrueSet.IsKey(i)) continue;
		for (int j = 0; j < GNodeNum; j++) {
			if (GTrueSet.IsKey(j)) continue;
			if (C[i][j] >= ColThH.GetDat(j)) continue;
			if (C[i][j] >= RowThH.GetDat(i)) continue;
			if (C[i][j] < ColMinThH.GetDat(j) || C[i][j] < RowMinThH.GetDat(i))
				C[i][j] = 0;
		}
	}

	int ErrorCnt = 0;
	// Normalize according to column
	for (int j = 0; j < GNodeNum; j++) {
		if (GTrueSet.IsKey(j)) continue;
		double Z = 0;
		for (int i = 0; i < FNodeNum; i++) {
			if (FTrueSet.IsKey(i)) continue;
			Z += C[i][j];
		}
		if (Z < 1e-5) ErrorCnt += 1;
		for (int i = 0; i < FNodeNum; i++)
			C[i][j] /= Z;
	}
	printf("ErrorCnt = %d\n", ErrorCnt);
		

	// Normalize according to row
	for (int i = 0; i < FNodeNum; i++) {
		if (FTrueSet.IsKey(i)) continue;
		double Z = 0;
		for (int j = 0; j < GNodeNum; j++) {
			if (GTrueSet.IsKey(j)) continue;
			Z += C[i][j];
		}
		for (int j = 0; j < GNodeNum; j++)
			C[i][j] /= Z;
	}
}

void SaveRData(TStr Pref, int Num) {
	TStr FNm = TStr::Fmt("%s-R%d.bin", Pref.CStr(), Num);
	FILE* FLog = fopen(FNm.CStr(), "w");
	for (int i = 0; i < FNodeNum; i++) {
		for (int j = 0; j < GNodeNum; j++)
			fprintf(FLog, "%lg\t", CurR[i][j]);
		fprintf(FLog, "\n");
	}
	fclose(FLog);
}

void SaveCData(TStr Pref) {
	TStr FNm = TStr::Fmt("%s-C.bin", Pref.CStr());
	FILE* FLog = fopen(FNm.CStr(), "w");
	for (int i = 0; i < FNodeNum; i++) {
		for (int j = 0; j < GNodeNum; j++)
			fprintf(FLog, "%lg\t", C[i][j]);
		fprintf(FLog, "\n");
	}
	fclose(FLog);
}

void LoadRData(TStr Pref, int Num) {
	TStr FNm = TStr::Fmt("%s-R%d.bin", Pref.CStr(), Num);
	FILE* FLog = fopen(FNm.CStr(), "r");
	for (int i = 0; i < FNodeNum; i++) {
		for (int j = 0; j < GNodeNum; j++)
			fscanf(FLog, "%lf", &CurR[i][j]);
	}	
	fclose(FLog);
}

void LoadCData(TStr Pref) {
	TStr FNm = TStr::Fmt("%s-C.bin", Pref.CStr());
	printf("Loading from %s\n", FNm.CStr());
	FILE* FLog = fopen(FNm.CStr(), "r");
	for (int i = 0; i < FNodeNum; i++) {
		for (int j = 0; j < GNodeNum; j++)
			fscanf(FLog, "%lf", &C[i][j]);
	}	
	double AvgDeg = 0, TotCnt = 0;
	for (int i = 0; i < FNodeNum; i++) {
		int Deg = 0;
		for (int j = 0; j < GNodeNum; j++) {
			if (C[i][j] > 1e-8)
				Deg += 1;
		}
		if (Deg < 10) continue;
		AvgDeg += Deg;
		TotCnt += 1;
	}
	AvgDeg /= TotCnt;
	printf("Avg Deg in C Matrix is %.4lf\n", AvgDeg);
	fclose(FLog);
}

void IsoRank(double alpha, int PIter, TStr Pref) {
	// Normalize C
	double TotC = 0;
	for (int i = 0; i < FNodeNum; i++)
		for (int j = 0; j < GNodeNum; j++)
			TotC += C[i][j];
	TotC /= FNodeNum;
	for (int i = 0; i < FNodeNum; i++)
		for (int j = 0; j < GNodeNum; j++)
			C[i][j] /= TotC;

	IAssert(PIter >= FirstIter);
	for (int Iter = 0; Iter < 12; Iter++) {
		if (PIter == FirstIter && Iter < SecondIter) continue;
		if (PIter == FirstIter && Iter == SecondIter) {
			LoadRData(Pref, Iter);
		} else {
			SaveRData(Pref, Iter);
		}
		printf("Start IsoRank Iteration %d\n", Iter);
		// NextR = A * CurR
		int OmpCnt = 0;
		#pragma omp parallel for schedule(dynamic, 10)
		for (int i = 0; i < FNodeNum; i++) {
			if (OmpCnt++ % 100 == 0) printf("i = %d\t", OmpCnt);
			for (int j = 0; j < GNodeNum; j++) {
				double sum = 0;
				for (int Nbi = 0; Nbi < LenF[i]; Nbi++)
					for (int Nbj = 0; Nbj < LenF[j]; Nbj++)
						sum += CurR[NbF[i][Nbi]][NbG[j][Nbj]] * NbSF[i][Nbi] * NbSG[j][Nbj];
				NextR[i][j] = sum;
			}
		}

		// NextR = alpha * NextR + (1 - alpha) * C
		double TotR = 0;
		for (int i = 0; i < FNodeNum; i++) 
			for (int j = 0; j < GNodeNum; j++)
				TotR += NextR[i][j];
		TotR /= FNodeNum;
		for (int i = 0; i < FNodeNum; i++) 
			for (int j = 0; j < GNodeNum; j++)
				NextR[i][j] = alpha * NextR[i][j]/TotR  + (1-alpha) * C[i][j];
		double RDif = ComputeRelativeRDif();

		// CurR = NextR
		for (int i = 0; i < FNodeNum; i++)
			for (int j = 0; j < GNodeNum; j++)
				CurR[i][j] = NextR[i][j] / DegF[i] / DegG[j];
		if (RDif < 1e-4) break;
	}
}

void MatchGraph(PNGraph& F, PNGraph& G, PNGraph& M, PNGraph& RetM) {
	printf("Start matching graphs\n");
	ComputeFDegNb(F, "LOG_FDegNb.txt");
	ComputeGDegNb(G, "LOG_GDegNb.txt");
//	DumpToMatlab(F, G, M);
	double alpha = 0.5;
	InitializeC(M);
	InitializeR();
	double thresh = 0.4;
	for (int Iter = 0; Iter < 9; Iter++) {
		TStr Pref = TStr::Fmt("Data%d", Iter);
		if (Iter < FirstIter) {thresh = thresh * 0.4; continue;}
		if (Iter == FirstIter && SecondIter >= 0)
			LoadCData(Pref);
		else
			SaveCData(Pref);
		printf("-------------------------- Iteration %d ------------------------------\n", Iter);

		IsoRank(alpha, Iter, Pref);
		ComputeC(M, thresh, Pref);
		thresh = thresh * 0.4;
	}
}

int main() {
	TSecTm BegTm = TSecTm::GetCurTm();
	PNGraph F = PNGraph::New();
	PNGraph G = PNGraph::New();
	PNGraph M = PNGraph::New();
	LoadGraph(F, "F.txt");
	LoadGraph(G, "G.txt");
	LoadGraph(M, "M.txt");
	for (TNGraph::TNodeI NI = M->BegNI(); NI < M->EndNI(); NI++) {
		if (NI.GetInDeg() == 0) {
			FSet.AddKey(NI.GetId());
			FMatchH.AddDat(NI.GetId()) = NI.GetOutNId(0);
		}	else {
			GSet.AddKey(NI.GetId());
			GMatchH.AddDat(NI.GetId()) = NI.GetInNId(0);
		}
	}
	//TESTDumpTrueSetInfo(F, G, M);
	GetMatchNum(F, G, M);
	PNGraph RetM = PNGraph::New();
	MatchGraph(F, G, M, RetM);
	GetMatchNum(F, G, RetM);

  TSecTm EndTm = TSecTm::GetCurTm();
  double usedTime = EndTm.GetAbsSecs() - BegTm.GetAbsSecs();
  printf("Total execution time : %02dh%02dm%02ds\n", int(usedTime)/3600, (int(usedTime)%3600)/60, int(usedTime)%60);
	return 0;
}


void TESTDumpEdgeMatchInfo(PNGraph& F, PNGraph& G, TIntH& FMatchH, TIntH& GMatchH) {
	int TotTrueCnt = 0, FwdEdgeCnt = 0, BackEdgeCnt = 0, OverlapEdgeCnt = 0;
	for (TNGraph::TEdgeI EI = F->BegEI(); EI < F->EndEI(); EI++) {
		if (!FMatchH.IsKey(EI.GetSrcNId()) || !FMatchH.IsKey(EI.GetDstNId())) continue;
		int FSrcId = EI.GetSrcNId(), FDstId = EI.GetDstNId();
		int GSrcId = FMatchH.GetDat(FSrcId), GDstId = FMatchH.GetDat(FDstId);
		TotTrueCnt += 1;
		if (G->IsEdge(GSrcId, GDstId))
			FwdEdgeCnt += 1;
		if (G->IsEdge(GDstId, GSrcId))
			BackEdgeCnt += 1;
		if (G->IsEdge(GSrcId, GDstId) && G->IsEdge(GDstId, GSrcId))
			OverlapEdgeCnt += 1;
	}
	printf("Total Edge = %d, Fwd = %d, Back = %d, Overlap = %d\n", TotTrueCnt, FwdEdgeCnt, BackEdgeCnt, OverlapEdgeCnt);
}

void TESTDumpTrueSetInfo(PNGraph& F, PNGraph& G, PNGraph& M) {
	ComputePageRank(F, FPageRankH);
	ComputePageRank(G, GPageRankH);
	FILE* FTMP = fopen("TEMP_log.txt", "w");
	for (int i = 0; i < FSet.Len(); i++) {
		int FId = FSet[i];
		int GId = FMatchH.GetDat(FId);
		IAssert(GMatchH.GetDat(GId) == FId);
		fprintf(FTMP, "F : Id = %d, InDeg = %d, OutDeg = %d, PageRank = %f\n", FId, F->GetNI(FId).GetInDeg(), F->GetNI(FId).GetOutDeg(), FPageRankH.AddDat(FId).Val);
		fprintf(FTMP, "G : Id = %d, InDeg = %d, OutDeg = %d, PageRank = %f\n\n", GId, G->GetNI(GId).GetInDeg(), G->GetNI(GId).GetOutDeg(), GPageRankH.AddDat(GId).Val);
	}
	fclose(FTMP);
	TESTDumpEdgeMatchInfo(F, G, FMatchH, GMatchH);
	TESTDumpEdgeMatchInfo(G, F, GMatchH, FMatchH);
}

void GetMatchNum(PNGraph& F, PNGraph& G, PNGraph& M) {
	int TotNum = 0, FBackNum = 0, GBackNum = 0;
	for (TNGraph::TEdgeI EI = F->BegEI(); EI < F->EndEI(); EI++) {
		int FSrcId = EI.GetSrcNId(), FDstId = EI.GetDstNId();
		IAssert(FSrcId != FDstId);
		if (!M->IsNode(FSrcId)) continue;
		if (!M->IsNode(FDstId)) continue;
		int GSrcId = M->GetNI(FSrcId).GetOutNId(0);
		int GDstId = M->GetNI(FDstId).GetOutNId(0);
		if (G->IsEdge(GSrcId, GDstId))
			TotNum += 1;
		else if (!F->IsEdge(FDstId, FSrcId) && G->IsEdge(GDstId, GSrcId))
			FBackNum += 1;
	}

	for (TNGraph::TEdgeI EI = G->BegEI(); EI < G->EndEI(); EI++) {
		int GSrcId = EI.GetSrcNId(), GDstId = EI.GetDstNId();
		IAssert(GSrcId != GDstId);
		if (!M->IsNode(GSrcId)) continue;
		if (!M->IsNode(GDstId)) continue;
		int FSrcId = M->GetNI(GSrcId).GetInNId(0);
		int FDstId = M->GetNI(GDstId).GetInNId(0);
		if (!F->IsEdge(FSrcId, FDstId) && !G->IsEdge(GDstId, GSrcId) && F->IsEdge(FDstId, FSrcId))
			GBackNum += 1;
	}
	printf("\n\nSTATS:\nTotNum = %d, FBackNum = %d GBackNum = %d\n\n", TotNum, FBackNum, GBackNum);
}

void LoadGraph(PNGraph& G, TStr FNm) {
	TFIn FIn(FNm);
	for (TChA Ln; FIn.GetNextLn(Ln);) {
		if (Ln.IsPrefix("#")) continue;
		if (Ln.Len() <= 1) continue;
		int idx = 0;
		while (Ln[idx] != '\t') idx++;
		Ln[idx] = 0;
		int SrcId, DstId;
		SrcId = atoi(Ln.CStr());
		DstId = atoi(Ln.CStr() + idx + 1);
		if (!G->IsNode(SrcId))
			G->AddNode(SrcId);
		if (!G->IsNode(DstId))
			G->AddNode(DstId);
		G->AddEdge(SrcId, DstId);
	}
	TSnap::PrintInfo(G);
}

void ComputeAvgDist(const PNGraph& G, TIntFltH& AvgDistH, TIntFltH& FwdNumH, TIntFltH& BackNumH, const TIntSet& IdSet) {
	printf("Computing average distance\n");
	TIntV IdV;
	for (TNGraph::TNodeI NI = G->BegNI(); NI < G->EndNI(); NI++)
		IdV.Add(NI.GetId());
	TFltV AvgDistRecV(IdV.Len());

	#pragma omp parallel for
	for (int i = 0; i < IdV.Len(); i++) {
		int Id = IdV[i];
		TIntV Queue; 
		TIntH DistH; DistH.AddDat(Id) = 0;
		int idx = -1;
		Queue.Add(Id); 
		while (idx < Queue.Len() - 1) {
			idx += 1;
			int CurId = Queue[idx];
			TNGraph::TNodeI CurNI = G->GetNI(CurId);
			for (int e = 0; e < CurNI.GetInDeg(); e++) {
				if (DistH.IsKey(CurNI.GetInNId(e))) continue;
				int NextId = CurNI.GetInNId(e);
				DistH.AddDat(NextId) = DistH.GetDat(CurId) + 1;
				Queue.Add(NextId);
			}
			for (int e = 0; e < CurNI.GetOutDeg(); e++) {
				if (DistH.IsKey(CurNI.GetOutNId(e))) continue;
				int NextId = CurNI.GetOutNId(e);
				DistH.AddDat(NextId) = DistH.GetDat(CurId) + 1;
				Queue.Add(NextId);
			}
		}
		double avgdist = 0;
		for (int j = 0; j < IdSet.Len(); j++) 
			avgdist += DistH.GetDat(IdSet[j]).Val;
		avgdist /= IdSet.Len();
		AvgDistRecV[i] = avgdist;
	}
	for (int i = 0; i < IdV.Len(); i++) {
		int Id = IdV[i]; double avgdist = AvgDistRecV[i];
		AvgDistH.AddDat(Id) = avgdist;
	}

	TFltV FwdNumRecV(IdV.Len());
	// Forward Number
	printf("Computing forward number\n");
	#pragma omp parallel for
	for (int i = 0; i < IdV.Len(); i++) {
		int Id = IdV[i];
		TIntV Queue; 
		TIntH DistH; DistH.AddDat(Id) = 0;
		int idx = -1;
		Queue.Add(Id); 
		while (idx < Queue.Len() - 1) {
			idx += 1;
			int CurId = Queue[idx];
			if (DistH.GetDat(CurId) > 1) break;
			TNGraph::TNodeI CurNI = G->GetNI(CurId);
			for (int e = 0; e < CurNI.GetOutDeg(); e++) {
				if (DistH.IsKey(CurNI.GetOutNId(e))) continue;
				int NextId = CurNI.GetOutNId(e);
				DistH.AddDat(NextId) = DistH.GetDat(CurId) + 1;
				Queue.Add(NextId);
			}
		}
		double fwdnum = 0;
		for (int j = 0; j < idx; j++)
			if (IdSet.IsKey(Queue[j])) fwdnum += 1;
		FwdNumRecV[i] = fwdnum;
	}
	for (int i = 0; i < IdV.Len(); i++) {
		int Id = IdV[i];
		double fwdnum = FwdNumRecV[i];
		FwdNumH.AddDat(Id) = fwdnum;
	}

	TFltV BackNumRecV(IdV.Len());
	// Backward Number
	printf("Computing backward number\n");
	#pragma omp parallel for
	for (int i = 0; i < IdV.Len(); i++) {
		int Id = IdV[i];
		TIntV Queue; 
		TIntH DistH; DistH.AddDat(Id) = 0;
		int idx = -1;
		Queue.Add(Id); 
		while (idx < Queue.Len() - 1) {
			idx += 1;
			int CurId = Queue[idx];
			if (DistH.GetDat(CurId) > 1) break;
			TNGraph::TNodeI CurNI = G->GetNI(CurId);
			for (int e = 0; e < CurNI.GetInDeg(); e++) {
				if (DistH.IsKey(CurNI.GetInNId(e))) continue;
				int NextId = CurNI.GetInNId(e);
				DistH.AddDat(NextId) = DistH.GetDat(CurId) + 1;
				Queue.Add(NextId);
			}
		}
		double backnum = 0;
		for (int j = 0; j < idx; j++)
			if (IdSet.IsKey(Queue[j])) backnum += 1;
		BackNumRecV[i] = backnum;
	}
	for (int i = 0; i < IdV.Len(); i++) {
		int Id = IdV[i];
		double backnum = BackNumRecV[i];
		BackNumH.AddDat(Id) = backnum;
	}
}

void PolishGraph(PNGraph& GG, const PNGraph& G) {
	for (TNGraph::TNodeI NI = G->BegNI(); NI < G->EndNI(); NI++) {
		GG->AddNode(NI.GetId());
	}
	for (TNGraph::TEdgeI EI = G->BegEI(); EI < G->EndEI(); EI++) {
		GG->AddEdge(EI.GetSrcNId(), EI.GetDstNId());
	}

	while (1) {
		TIntV DelNodeV;
		for (TNGraph::TNodeI NI = GG->BegNI(); NI < GG->EndNI(); NI++) 
			if (NI.GetInDeg() == 0 || NI.GetOutDeg() == 0)
				DelNodeV.Add(NI.GetId());
		if (DelNodeV.Len() == 0) break;
		for (int i = 0; i < DelNodeV.Len(); i++)
			GG->DelNode(DelNodeV[i]);
	}
	TSnap::PrintInfo(GG);
}

void ComputePageRank(const PNGraph& GG, TIntFltH& PageRankH) {
	printf("Computing page rank scores\n");

	PNGraph G = PNGraph::New();
	PolishGraph(G, GG);

	TIntV IdV;
	for (TNGraph::TNodeI NI = G->BegNI(); NI < G->EndNI(); NI++)
		IdV.Add(NI.GetId());

	TIntFltH CurH;
	for (int i = 0; i < IdV.Len(); i++)
		CurH.AddDat(IdV[i]) = 1.0; 
	
	int IterNum = 100;
	for (int Iter = 0; Iter < IterNum; Iter++) {
		if (Iter % 50 == 0) printf("%d out of %d completed\n", Iter, IterNum);
		TFltV NextV(IdV.Len());
		#pragma omp parallel for
		for (int i = 0; i < IdV.Len(); i++) {
			int Id = IdV[i];
			TNGraph::TNodeI NI = G->GetNI(Id);
			for (int e = 0; e < NI.GetInDeg(); e++) {
				TNGraph::TNodeI NeighborNI = G->GetNI(NI.GetInNId(e));
				NextV[i] += CurH.GetDat(NeighborNI.GetId()).Val / 1.0 / NeighborNI.GetOutDeg();
			}
		}		
		for (int i = 0; i < IdV.Len(); i++)
			CurH.GetDat(IdV[i]) = NextV[i];
	}
	PageRankH = CurH;
}

void DumpToMatlab(PNGraph& F, PNGraph& G, PNGraph& M) {
	printf("Dumping info for MATLAB\n");
	FILE* FMAT = fopen("FMAT.txt", "w");
	fprintf(FMAT, "%d\t%d\n", FNodeNum, GNodeNum);
	int EdgeNum = 0;
	for (int i = 0; i < FNodeNum; i++) {
		for (int j = 0; j < FNodeNum; j++) {
			if (F->IsEdge(FNodeIdH.GetKey(i), FNodeIdH.GetKey(j))) {
				fprintf(FMAT, "1\t");
				EdgeNum += 1;
			}	else {
				fprintf(FMAT, "0\t");
			}
		}
		fprintf(FMAT, "\n");
	}
	printf("EdgeNum = %d\n", EdgeNum);
	fclose(FMAT);

	FILE* GMAT = fopen("GMAT.txt", "w");
	EdgeNum = 0;
	fprintf(GMAT, "%d\t%d\n", FNodeNum, GNodeNum);
	for (int i = 0; i < GNodeNum; i++) {
		for (int j = 0; j < GNodeNum; j++) {
			if (G->IsEdge(GNodeIdH.GetKey(i), GNodeIdH.GetKey(j))) {
				fprintf(GMAT, "1\t");
				EdgeNum += 1;
			} else {
				fprintf(GMAT, "0\t");
			}
		}
		fprintf(GMAT, "\n");
	}
	printf("EdgeNum = %d\n", EdgeNum);
	fclose(GMAT);

	FILE* MMAT = fopen("MMAT.txt", "w");
	fprintf(MMAT, "%d\n", 800);
	for (TNGraph::TEdgeI EI = M->BegEI(); EI < M->EndEI(); EI++) {
		int SrcN = FNodeIdH.GetDat(EI.GetSrcNId());
		int DstN = GNodeIdH.GetDat(EI.GetDstNId());
		fprintf(MMAT, "%d\t%d\n", SrcN, DstN);
	}
	fclose(MMAT);
	printf("DONE\n");
}

